data {
  int<lower=0> P;
  int<lower=0> N;
  int<lower=0> I;
  int<lower=0> C;
  int<lower=0> E;
  
  int<lower=0> F1;
  int<lower=0> F2;
  
  int<lower=1,upper=11> p[N];
  
  int<lower=0,upper=C> c_id[N];
  int<lower=0,upper=E> e_id[N];
  int<lower=0,upper=P> p_id[N];
  int<lower=0,upper=I> i_id[N];
  
  int<lower=0,upper=I> f1[I-F1];
  int<lower=0,upper=I> f2[I-F2];
  
  int<lower=0,upper=I> f1z[1,F1];
  int<lower=0,upper=I> f2z[1,F2];
  
  int<lower=0,upper=1> free_experts;
  int<lower=0,upper=1> party_issue_shocks;
  int<lower=0,upper=1> country_issue_shocks;
  int<lower=0,upper=1> country_shocks;
}

parameters {
  matrix[P,2] theta_free;
  
  ordered[10] mu_alpha[I];
  simplex[I] scales_alpha;
  cholesky_factor_corr[I] L_alpha;
  matrix[I,C] res_alpha_r;
  
  vector[I-F1] mu_beta_1;
  vector[I-F2] mu_beta_2;
  simplex[I-F1] scales_beta_1;
  simplex[I-F2] scales_beta_2;
  cholesky_factor_corr[I-F1] L_beta_1;
  cholesky_factor_corr[I-F2] L_beta_2;
  matrix[I-F1,C] res_beta_r1;
  matrix[I-F2,C] res_beta_r2;
  
  simplex[C] scale_beta_free;
  vector[C] shift_alpha_free;
  real<lower=0> hyp_p[6];

  simplex[E] exp_scale_free;
  matrix[I,P] idio_r;
}

transformed parameters {
  matrix[I-F1,C] res_beta_1;
  matrix[I-F2,C] res_beta_2;
  matrix[I,C] beta_1;
  matrix[I,C] beta_2;
  
  matrix[I,C] res_alpha;
  
  vector<lower=0>[E] exp_scale;
  
  matrix[I,P] idio;
  matrix[P,2] theta;
  
  vector<lower=0>[C] scale_beta;
  vector[C] shift_alpha;
  
  vector[N] ystar;
  
  for (d in 1:2){
    theta[,d] = (theta_free[,d] - mean(to_vector(theta_free[,d]))) / sd(to_vector(theta_free[,d]));
  }
  
  {
    matrix[I-F1,C] res_beta_temp1 = quad_form_diag(L_beta_1 * L_beta_1', sqrt(1.0/hyp_p[1] * (scales_beta_1))) * res_beta_r1;
    matrix[I-F2,C] res_beta_temp2 = quad_form_diag(L_beta_2 * L_beta_2', sqrt(1.0/hyp_p[1] * (scales_beta_2))) * res_beta_r2;
    matrix[I,C] res_alpha_temp = quad_form_diag(L_alpha * L_alpha', sqrt(1.0/hyp_p[2] * (scales_alpha))) * res_alpha_r;
    
    for (i in 1:(I-F1)) res_beta_1[i,] = res_beta_temp1[i,] - mean(res_beta_temp1[i,]);
    for (i in 1:(I-F2)) res_beta_2[i,] = res_beta_temp2[i,] - mean(res_beta_temp2[i,]);
    if (country_issue_shocks == 1){
      for (i in 1:(I)) res_alpha[i,] = res_alpha_temp[i,] - mean(res_alpha_temp[i,]);
    } else {
      for (i in 1:(I)) res_alpha[i,] = to_row_vector(rep_vector(0.0, C));
    }
    
    for (i in 1:(I-F1)){
      for (c in 1:C){
        if (country_issue_shocks == 1){
          beta_1[f1[i],c] = mu_beta_1[i] + res_beta_1[i,c];
        } else {
          beta_1[f1[i],c] = mu_beta_1[i];
        }
      }
    }
    
    for (i in 1:(I-F2)){
      for (c in 1:C){
        if (country_issue_shocks == 1){
          beta_2[f2[i],c] = mu_beta_2[i] + res_beta_2[i,c];
        } else {
          beta_2[f2[i],c] = mu_beta_2[i];
        }
      }
    }
    
    for (i in 1:(F1)) beta_1[f1z[1,i],] = rep_row_vector(0.0, C);
    for (i in 1:(F2)) beta_2[f2z[1,i],] = rep_row_vector(0.0, C);
  }
  
  if (free_experts == 0){
    exp_scale = rep_vector(1.0, E);
  } else {
    exp_scale = E * exp_scale_free;
  }
  
  if (country_shocks == 0){
    scale_beta = rep_vector(1.0, C);
    shift_alpha = rep_vector(0.0, C);
  } else {
    scale_beta = C * scale_beta_free;
    shift_alpha = (hyp_p[6] * shift_alpha_free) - mean(hyp_p[6] * shift_alpha_free);
  }
  
  if (party_issue_shocks == 0){
    for (i in 1:I) idio[i,] = rep_row_vector(0.0, P);
  } else {
    for (i in 1:I) idio[i,] = hyp_p[5] * idio_r[i,] - mean(hyp_p[5] * idio_r[i,]);
  }
  
  for (n in 1:N){
    ystar[n] = (scale_beta[c_id[n]] * (theta[p_id[n],1] * beta_1[i_id[n],c_id[n]] + theta[p_id[n],2] * beta_2[i_id[n],c_id[n]] + idio[i_id[n],p_id[n]])) / exp_scale[e_id[n]];
  }
}

model {
  vector[N] log_lik;
  for (i in 1:I) (mu_alpha[i]) ~ std_normal();
  
  mu_beta_1 ~ std_normal();
  mu_beta_2 ~ std_normal();
  
  to_vector(res_beta_r1) ~ std_normal();
  to_vector(res_beta_r2) ~ std_normal();
  to_vector(res_alpha_r) ~ std_normal();
  
  scales_beta_1 ~ dirichlet(rep_vector(1, I-F1));
  scales_beta_2 ~ dirichlet(rep_vector(1, I-F2));
  scales_alpha ~ dirichlet(rep_vector(1, I));
  
  L_beta_1 ~ lkj_corr_cholesky(4);
  L_beta_2 ~ lkj_corr_cholesky(4);
  L_alpha ~ lkj_corr_cholesky(4);
    
  to_vector(idio_r) ~ std_normal();
  
  to_vector(theta_free) ~ std_normal();

  hyp_p ~ std_normal();

  exp_scale_free ~ dirichlet(rep_vector(1/hyp_p[3], E));
  
  scale_beta_free ~ dirichlet(rep_vector(1/hyp_p[4], C));
  shift_alpha_free ~ std_normal();
  
  for (n in 1:N){
    if (i_id[n] < 17) {
      log_lik[n] = ordered_probit_lupmf(p[n] | ystar[n],
      (mu_alpha[i_id[n]] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]]);
    } else {
      log_lik[n] = ordered_probit_lupmf(p[n] | ystar[n],
      (mu_alpha[i_id[n]][1:6] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]]);
    }
  }
  target += sum(log_lik);
}
generated quantities{
  corr_matrix[I-F1] LL_beta1 = L_beta_1*L_beta_1';
  corr_matrix[I-F2] LL_beta2 = L_beta_2*L_beta_2';
  corr_matrix[I] LL_alpha = L_alpha*L_alpha';
  vector[N] log_lik;
  real pred = 0.0;
  real p_share; // share of correct votes
  for (n in 1:N){
    if (i_id[n] < 17) {
      log_lik[n] = ordered_probit_lpmf(p[n] | ystar[n],
      (mu_alpha[i_id[n]] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]]);
      if (p[n] == ordered_probit_rng(ystar[n], (mu_alpha[i_id[n]] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]])) {
        pred += 1.00;
        }
    } else {
      log_lik[n] = ordered_probit_lpmf(p[n] | ystar[n],
      (mu_alpha[i_id[n]][1:6] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]]);
      if (p[n] == ordered_probit_rng(ystar[n], (mu_alpha[i_id[n]] + shift_alpha[c_id[n]] + res_alpha[i_id[n],c_id[n]]) / exp_scale[e_id[n]])) {
        pred += 1.00;
        }
        }
      } 
    
  p_share = pred / (N*1.00);
}
